#include "matrix.h"

//
// Allocate a matrix of size (m, n)
//
double** allocate_matrix(int m, int n)
{
	double** M = new double*[m];
	for (int i = 0; i < m; i++)
		M[i] = new double[n];
	return M;
}

//
// Free a matrix of size (m, n)
//
void free_matrix(double** M, int m, int n)
{
	for (int i = 0; i < m; i++)
		delete[] M[i];
	delete[] M;
}

//
// Copy a matrix of size (m, n) to another matrix
//
double** copy_matrix(double** M, int m, int n)
{
	double** C = allocate_matrix(m, n);
	for (int i = 0; i < m; i++)
		for (int j = 0; j < n; j++)
			C[i][j] = M[i][j];
	return C;
}

//
// Returns the minor of the element (r, c) in the matrix M
//
void minor(double** M, int m, int n, int r, int c, double** minor)
{
	int rr = 0;
	int cc = 0;
	for (int i = 0; i < m; i++)
	{
		if (i != r)
		{
			for (int j = 0; j < n; j++)
				if (j != c)
					minor[rr][cc++] = M[i][j];
			rr++;
			cc = 0;
		}
	}
}

//
// Compute the determinent of the matrix
//
double determinent(double** M, int m, int n)
{
	if (m != n)
		return 0;

	if (m == 2)
		return M[0][0] * M[1][1] - M[0][1] * M[1][0];
	else
	{
		double det_all = 0.0;
		for (int c = 0; c < n; c++)
		{
			double** m_c = allocate_matrix(m - 1, n - 1);
			minor(M, m, n, 0, c, m_c);
			double det = determinent(m_c, m - 1, n - 1);
			det *= M[0][c];
			det *= (c % 2 == 0) ? 1 : -1;
			det_all += det;
			free_matrix(m_c, m - 1, n - 1);
		}
		return det_all;
	}

	return 0.0;
}

// 
// Solve the linear system using the cramers rule
//
void solve_cramer(double** A, int m, int n, double** B, double** X)
{
	double det_A = determinent(A, m, n);
	for (int c = 0; c < n; c++)
	{
		double** A_c = copy_matrix(A, m, n);
		for (int r = 0; r < m; r++)
			A_c[r][c] = B[r][0];
		double det_A_c = determinent(A_c, m, n);
		X[c][0] = det_A_c / det_A;
		free_matrix(A_c, m, n);
	}
}